Multispecies survey specifications

Species surveyed

Our initial species selection includes 11 single species groups from the NOBA model. We looked at diet overlap here. Below we will double check survey census outputs against the true values, and then test some survey specifications that may be reasonable for multispecies model testing. Finally we will apply survey functions to the detailed diet data.

fgs <- load_fgs(here("atlantisoutput", "NOBA_March_2020"), "nordic_groups_v04.csv")

lname <- data.frame(Latin = c("*Hippoglossoides platessoides*",
                              "*Reinhardtius hippoglossoides*",
                              "*Scomber scombrus*",
                              "*Melongrammus aeglefinus*",
                              "*Pollachius virens*",
                              "*Sebastes mentella*",
                              "*Micromesistius poutassou*",
                              "*Clupea harengus*",
                              "*Gadus morhua*",
                              "*Boreogadus saida*",
                              "*Mallotus villosus*"),
                    Code = c("LRD", "GRH", "MAC", "HAD", "SAI", "RED", 
                             "BWH", "SSH", "NCO", "PCO", "CAP")
)

sppsubset <- merge(fgs, lname, all.y = TRUE)
spptable <- sppsubset %>% 
  arrange(Index) %>%
  select(Name, Long.Name, Latin)

knitr::kable(spptable, col.names = c("Model name", "Full name", "Latin name"))

Model name	Full name	Latin name
Long_rough_dab	Long rough dab	Hippoglossoides platessoides
Green_halibut	Greenland halibut	Reinhardtius hippoglossoides
Mackerel	Mackerel	Scomber scombrus
Haddock	Haddock	Melongrammus aeglefinus
Saithe	Saithe	Pollachius virens
Redfish	Redfish	Sebastes mentella
Blue_whiting	Blue whiting	Micromesistius poutassou
Norwegian_ssh	Norwegian spring spawning herring	Clupea harengus
North_atl_cod	Northeast Atlantic cod	Gadus morhua
Polar_cod	Polar cod	Boreogadus saida
Capelin	Capelin	Mallotus villosus

Generating the census dataset (October 2020)

Note that survey config files must now name the survey.

Unchanging config files

NOBA2config.R looks like this (adjusted from Alfonso’s original):

d.name <- here("atlantisoutput","NOBA_march_2020")
functional.groups.file <- "nordic_groups_v04.csv"
biomass.pools.file <- "nordic_biol_v23.nc"
biol.prm.file <- "nordic_biol_incl_harv_v_007_3.prm"
box.file <- "Nordic02.bgm"
initial.conditions.file <- "nordic_biol_v23.nc"
run.prm.file <- "nordic_run_v01.xml"
scenario.name <- "nordic_runresults_01"
bioind.file <- "nordic_runresults_01BiomIndx.txt"
catch.file <- "nordic_runresults_01Catch.txt"
annage <- TRUE
fisheries.file <- "NoBAFisheries.csv"

omdimensions.R standardizes timesteps, etc. (this is part of atlantisom and should not need to be changed by the user):

#survey species inherited from omlist_ss
survspp <- omlist_ss$species_ss
# survey season and other time dimensioning parameters
# generalized timesteps all models
noutsteps <- omlist_ss$runpar$tstop/omlist_ss$runpar$outputstep
timeall <- c(0:noutsteps)
stepperyr <- if(omlist_ss$runpar$outputstepunit=="days") 365/omlist_ss$runpar$toutinc
midptyr <- round(median(seq(0,stepperyr)))

# model areas, subset in surveyconfig
allboxes <- c(0:(omlist_ss$boxpars$nbox - 1))

# fishery output: learned the hard way this can be different from ecosystem outputs
fstepperyr <- if(omlist_ss$runpar$outputstepunit=="days") 365/omlist_ss$runpar$toutfinc

# survey selectivity (agecl based)
sp_age <- omlist_ss$funct.group_ss[, c("Name", "NumCohorts", "NumAgeClassSize")]

# should return all age classes fully sampled (Atlantis output is 10 age groups per spp)
n_age_classes <- sp_age$NumCohorts
# changed below for multiple species NOTE survspp alphabetical; NOT in order of fgs!!
# this gives correct names
age_classes <- sapply(n_age_classes, seq)
names(age_classes)<-sp_age$Name

n_annages <- sp_age$NumCohorts * sp_age$NumAgeClassSize
# changed below for multiple species
annages <- sapply(n_annages, seq)
names(annages)<-sp_age$Name

Change these survey and fishery config files

mssurvey_spring.R and mssurvey_fall.R configure the fishery independent surveys (in this census test, surveys sample all model polygons in all years and have efficiency of 1 for all species, with no size selectivity):

# Default survey configuration here has a range of efficiencies and selectivities
# To emulate a range of species in a single multispecies survey
# Also now happens in "spring" and "fall"
# Need to define survey season, area, efficiency, selectivity

# Survey name
survey.name="BTS_spring_allbox_effic1"

#Atlantis model timestep corresponding to the true output--now from census_spec.R
timestep <- stepperyr #5

#Which atlantis timestep does the survey run in?--now from census_spec.R
# with 5 output steps per year, 0 is Jan-Feb-midMar, 1 is midMar-Apr-May, 
# 2 is June-July-midAug, 3 is midAug-Sept-Oct, 4 is Nov-Dec (ish)

survey_sample_time <- 1 # spring survey

#The last timestep to sample
total_sample <- noutsteps-1 #495

#Vector of indices of survey times to pull
survey_sample_full <- seq(survey_sample_time,
                          total_sample, by=timestep)

survtime <- survey_sample_full

# survey area
# should return all model areas
survboxes <- allboxes

# survey efficiency (q)
# should return a perfectly efficient survey 
surveffic <- data.frame(species=survspp,
                     efficiency=rep(1.0,length(survspp)))

# survey selectivity (agecl based)
# this is by age class, need to change to use with ANNAGEBIO output
#survselex <- data.frame(species=rep(names(age_classes), each=n_age_classes),
#                     agecl=rep(c(1:n_age_classes),length(survspp)),
#                     selex=rep(1.0,length(survspp)*n_age_classes))

# for annage output uses names(annages) NOT alphabetical survspp
survselex <- data.frame(species=rep(names(annages), n_annages), #  
                        agecl=unlist(sapply(n_annages,seq)),
                        selex=rep(1.0,sum(n_annages)))

# effective sample size needed for sample_fish
# this effective N is high but not equal to total for numerous groups
surveffN <- data.frame(species=survspp, effN=rep(100000, length(survspp)))

# survey index cv needed for sample_survey_xxx
# cv = 0.1
surv_cv <- data.frame(species=survspp, cv=rep(0.1,length(survspp)))

# length at age cv for input into calc_age2length function
# function designed to take one cv for all species, need to change to pass it a vector
lenage_cv <- 0.1

# max size bin for length estimation, function defaults to 150 cm if not supplied
maxbin <- 200

# Default survey configuration here has a range of efficiencies and selectivities
# To emulate a range of species in a single multispecies survey
# Also now happens in "spring" and "fall"
# Need to define survey season, area, efficiency, selectivity

# Survey name
survey.name="BTS_fall_allbox_effic1"

#Atlantis model timestep corresponding to the true output--now from census_spec.R
timestep <- stepperyr #5

#Which atlantis timestep does the survey run in?--now from census_spec.R
# with 5 output steps per year, 0 is Jan-Feb-midMar, 1 is midMar-Apr-May, 
# 2 is June-July-midAug, 3 is midAug-Sept-Oct, 4 is Nov-Dec (ish)

survey_sample_time <- 3 # fall survey

#The last timestep to sample
total_sample <- noutsteps-1 #495

#Vector of indices of survey times to pull
survey_sample_full <- seq(survey_sample_time,
                          total_sample, by=timestep)

survtime <- survey_sample_full

# survey area
# should return all model areas
survboxes <- allboxes

# survey efficiency (q)
# should return a perfectly efficient survey 
surveffic <- data.frame(species=survspp,
                     efficiency=rep(1.0,length(survspp)))

# survey selectivity (agecl based)
# this is by age class, need to change to use with ANNAGEBIO output
#survselex <- data.frame(species=rep(survspp, each=n_age_classes),
#                     agecl=rep(c(1:n_age_classes),length(survspp)),
#                     selex=rep(1.0,length(survspp)*n_age_classes))

# for annage output
survselex <- data.frame(species=rep(names(annages), n_annages), #  
                        agecl=unlist(sapply(n_annages,seq)),
                        selex=rep(1.0,sum(n_annages)))

# effective sample size needed for sample_fish
# this effective N is high but not equal to total for numerous groups
surveffN <- data.frame(species=survspp, effN=rep(100000, length(survspp)))

# survey index cv needed for sample_survey_xxx
# cv = 0.1
surv_cv <- data.frame(species=survspp, cv=rep(0.1,length(survspp)))

# length at age cv for input into calc_age2length function
# function designed to take one cv for all species, need to change to pass it a vector
lenage_cv <- 0.1

# max size bin for length estimation, function defaults to 150 cm if not supplied
maxbin <- 200

msfishery.R configures the fishery dependent data:

# Default fishery configuration here is a census
# Inherits species from input omlist_ss
fishspp <- omlist_ss$species_ss 

#Number of years of data to pull
nyears <- 50

#Atlantis initialization period in years
burnin <- 30

# fishery output: learned the hard way this can be different from ecosystem outputs
fstepperyr <- if(omlist_ss$runpar$outputstepunit=="days") 365/omlist_ss$runpar$toutfinc


# same time dimensioning parameters as in surveycensus.R
#Vector of indices of catch in numbers to pull (by timestep to sum)
fish_sample_full <- c(0:total_sample)  #total_sample defined in sardinesurvey.R
fish_burnin <- burnin*fstepperyr+1
fish_nyears <- nyears*fstepperyr
fish_times <- fish_sample_full[fish_burnin:(fish_burnin+fish_nyears-1)]
fish_timesteps <- seq(fish_times[fstepperyr], max(fish_times), by=fstepperyr) #last timestep
#fish_years <- unique(floor(fish_times/fstepperyr)+1) # my original
fish_years <- unique(floor(fish_times/fstepperyr)) #from Christine's new sardine_config.R

fishtime <- fish_times


# fishery sampling area
# should return all model areas, this assumes you see everything that it caught
fishboxes <- c(0:(omlist_ss$boxpars$nbox - 1))

# effective sample size needed for sample_fish
# this effective N is divided by the number of annual timesteps below, so 200 per time
# use as input to the length samples, ages can be a subset
fisheffN <- data.frame(species=survspp, effN=rep(1000, length(survspp)))

#this adjusts for subannual fishery output so original effN is for whole year
fisheffN$effN <- fisheffN$effN/fstepperyr 

# fishery catch cv can be used in sample_survey_biomass
# perfect observation
fish_cv <- data.frame(species=survspp, cv=rep(0.01,length(survspp)))

These survey and fishery configurations were used previously to generate saved census survey output.

Read in previously generated survey and fishery outputs and plot them against true outputs

Using read_savedsurvs() we get back the same objects generated by the wrapper functions om_index() and om_comps().

source(here("config", "NOBA2config.R"))
omlist_ss <- readRDS(file.path(d.name, paste0(scenario.name, "omlist_ss.rds")))
source(here("config/omdimensions.R"))

# add this paste0(d.name, "/mscensus1") and add option to name save directory in wrappers
survObsBiom <- read_savedsurvs(paste0(d.name, "/mscensus1"), 'survB')
age_comp_data <- read_savedsurvs(paste0(d.name, "/mscensus1"), 'survAge')
len_comp_data <- read_savedsurvs(paste0(d.name, "/mscensus1"), 'survLen')
wtage <- read_savedsurvs(paste0(d.name, "/mscensus1"), 'survWtage')
annage_comp_data <- read_savedsurvs(paste0(d.name, "/mscensus1"), 'survAnnAge')
annage_wtage <- read_savedsurvs(paste0(d.name, "/mscensus1"), 'survAnnWtage')

# no function needed for saved fishery output as long as there is only one fishery specification file

catchbio_ss <- readRDS(file.path(paste0(d.name, "/mscensus1"), paste0(scenario.name, "fishCatch.rds")))
catchbio_ss <- catchbio_ss[[1]] #in same time unit as surveys

fish_age_comp <- readRDS(file.path(paste0(d.name, "/mscensus1"), paste0(scenario.name, "fishObsAgeComp.rds")))
fish_age_comp <- fish_age_comp[[1]]

#add this to om_indices function so that this has years when read in
#fish_age_comp$time <- fish_age_comp$time/fstepperyr

fish_len_comp_data <- readRDS(file.path(paste0(d.name, "/mscensus1"), paste0(scenario.name, "fishObsLenComp.rds")))

#len_comp_data <- len_comp_data[[1]]
fish_len_comp_data <- fish_len_comp_data[[1]]

#add this to om_indices function so that this has years when read in
#fish_len_comp_data$time <- as.integer(floor(fish_len_comp_data$time/fstepperyr))

fishwtage <- readRDS(file.path(paste0(d.name, "/mscensus1"), paste0(scenario.name, "fishObsWtAtAge.rds")))
fishwtage <- fishwtage[[1]] #need to change time units?

fish_annage_comp <- readRDS(file.path(paste0(d.name, "/mscensus1"), paste0(scenario.name, "fishObsFullAgeComp.rds")))
fish_annage_comp <- fish_annage_comp[[1]]

#add this to om_indices function so that this has years when read in
#fish_annage_comp$time <- fish_annage_comp$time/fstepperyr

fish_annage_wtage <- readRDS(file.path(paste0(d.name, "/mscensus1"), paste0(scenario.name, "fishObsFullWtAtAge.rds")))
fish_annage_wtage <- fish_annage_wtage[[1]] #need to change time units?

Crunch a few things to get true age comps

# age classes
totNagecl <- omlist_ss$truenums_ss %>%
  group_by(species, agecl, time) %>%
  summarize(numAtAge = sum(atoutput))

totN <- totNagecl %>%
  group_by(species, time) %>%
  summarize(totN = sum(numAtAge))

trueNagecl <- merge(totNagecl, totN)

# annual ages
totNage <- omlist_ss$truenumsage_ss %>%
  group_by(species, agecl, time) %>%
  summarize(numAtAge = sum(atoutput))

#this should be exactly the same as totN but just making sure...
totN2 <- totNage %>%
  group_by(species, time) %>%
  summarize(totN = sum(numAtAge))

trueNage <- merge(totNage, totN2)

# fishery catch at age class
totCagecl <- omlist_ss$truecatchnum_ss %>%
  group_by(species, agecl, time) %>%
  summarize(numAtAge = sum(atoutput))

totCAAcl <- totCagecl %>%
  group_by(species, time) %>%
  summarize(totN = sum(numAtAge))

trueCAAcl <- merge(totCagecl, totCAAcl)

# fishery catch at annual age
totCage <- omlist_ss$truecatchage_ss %>%
  group_by(species, agecl, time) %>%
  summarize(numAtAge = sum(atoutput))

totCAA <- totCage %>%
  group_by(species, time) %>%
  summarize(totN = sum(numAtAge))

trueCAA <- merge(totCage, totCAA)

My plotting functions: too specialized to include in the package?

# plot biomass time series facet wrapped by species
plotB <- function(dat, truedat=NULL){
  
    ggplot() +
    geom_line(data=dat, aes(x=time/stepperyr,y=atoutput, color="Survey Biomass"), 
              alpha = 10/10) +
    {if(!is.null(truedat)) geom_line(data=truedat, aes(x=time/365,y=atoutput, color="True B"), alpha = 3/10)} + 
    theme_tufte() +
    theme(legend.position = "top") +
    xlab("model year") +
    ylab("tons") +
    labs(colour=scenario.name) +
    facet_wrap(~species, scales="free") 
  
}

# make a catch series function that can be split by fleet? this doesnt
# also note different time (days) from model timestep in all other output
plotC <- function(dat, truedat=NULL){
  
    ggplot() +
    geom_line(data=dat, aes(x=time/365,y=atoutput, color="Catch biomass"), 
              alpha = 10/10) +
    {if(!is.null(truedat)) geom_line(data=truedat, aes(x=time/365,y=atoutput, color="True Catch"), alpha = 3/10)} + 
    theme_tufte() +
    theme(legend.position = "top") +
    xlab("model year") +
    ylab("tons") +
    labs(colour=scenario.name) +
    facet_wrap(~species, scales="free") 
  
}

# plot length frequencies by timestep (one species)
plotlen <- function(dat, truedat=NULL){
  
    ggplot(dat, aes(upper.bins)) +
    geom_bar(aes(weight = atoutput)) +
    theme_tufte() +
    xlab("length (cm)") +
    ylab("number") +
    labs(subtitle = paste(scenario.name,
                          dat$species)) +
    facet_wrap(~time, ncol=6, scales="free_y")

}

# plot numbers at age by timestep (one species)
Natageplot <- function(dat, effN=1, truedat=NULL){
  ggplot() +
    geom_point(data=dat, aes(x=agecl, y=atoutput/effN, color="Est Comp")) +
    {if(!is.null(truedat)) geom_line(data=dat, aes(x=agecl, y=numAtAge/totN, color="True Comp"))} + 
    theme_tufte() +
    theme(legend.position = "bottom") +    
    xlab("age/agecl") +
    ylab("number") +
    labs(subtitle = paste(scenario.name,
                          dat$species)) + 
    facet_wrap(~time, ncol=6, scales="free_y")
}

# plot weight at age time series facet wrapped by species
wageplot <- function(dat, truedat=NULL){
  ggplot(dat, aes(time/stepperyr, atoutput)) +
    geom_line(aes(colour = factor(agecl))) +
    theme_tufte() +
    theme(legend.position = "bottom") +
    xlab("model year") +
    ylab("average individual weight (g)") +
    labs(subtitle = paste0(scenario.name)) +
    facet_wrap(c("species"), scales="free_y")
}
  
# compare N at age and C at age between standard and ANNAGE outputs
totNageplot <- function(dat, anndat){
  
    ggplot() +
    geom_line(data=dat, aes(x=time/stepperyr,y=totN, color="Tot N cohorts"), 
              alpha = 5/10) +
    geom_line(data=anndat, aes(x=time/stepperyr,y=totN, color="Tot N annage"), alpha = 5/10) + 
    theme_tufte() +
    theme(legend.position = "top") +
    xlab("model year") +
    ylab("N") +
    labs(colour=scenario.name) +
    facet_wrap(~species, scales="free") 
  
}

Visualize census surveys

We are using the full time series for biomass and weight at age, and a subsample of years 30-53 for compositions. In this test we see the effect of surveying in one season, with migratory species missed or only partially represented in some surveys, and with the overall surveyed biomass being at the lower or higher end of the range depending on the time of year. (Observation error is included in the survey biomass time series.) Samples for length and age assume a total of 100000 fish were randomly sampled for each species in each survey (surveffN defined in the survey config files).

Survey biomass

# compare with true output (all timesteps)
for(s in names(survObsBiom)){
  cat("  \n##### ",  s,"  \n")
  print(plotB(survObsBiom[[s]][[1]], omlist_ss$truetotbio_ss))
  cat("  \n")
}

BTS_fall_allbox_effic1

BTS_spring_allbox_effic1

# plots survey only
# for(s in names(survObsBiom)){
#   print(plotB(survObsBiom[[s]][[1]]))
# }

Survey length composition

for(s in names(len_comp_data)){
  cat("  \n##### ",  s,"  \n")
  lcompsub <- as.data.frame(len_comp_data[[s]][[1]]) %>% filter(time %in% c(150:270)) %>%
    group_by(species) %>%
    group_map(~ plotlen(.x), keep = TRUE)
  
  for(i in 1:length(lcompsub)) {
    print(lcompsub[[i]])
  }
  cat("  \n")
}

BTS_fall_allbox_effic1

BTS_spring_allbox_effic1

Survey age class composition (Atlantis age class)

for(s in names(age_comp_data)){
  cat("  \n##### ",  s,"  \n")
  acompsub <- as.data.frame(age_comp_data[[s]][[1]]) %>% filter(time %in% c(150:270)) %>%
    group_by(species) %>%
    left_join(., trueNagecl) %>%
    #group_map(~ Natageplot(.x), keep = TRUE) # plots only sampled age comp
    group_map(~ Natageplot(.x, effN = 100000, truedat = 1), keep = TRUE) # plots merged true age comp

  for(i in 1:length(acompsub)) {
    print(acompsub[[i]])
  }
  cat("  \n")
}

BTS_fall_allbox_effic1

BTS_spring_allbox_effic1

Survey age composition (annual ages)

for(s in names(annage_comp_data)){
  cat("  \n##### ",  s,"  \n")
  acompsub <- as.data.frame(annage_comp_data[[s]][[1]]) %>% filter(time %in% c(150:270)) %>%
    group_by(species) %>%
    left_join(., trueNage) %>%
    #group_map(~ Natageplot(.x), keep = TRUE) # plots only sampled age comp
    group_map(~ Natageplot(.x, effN = 100000, truedat = 1), keep = TRUE) # plots merged true age comp
  
  for(i in 1:length(acompsub)) {
    print(acompsub[[i]])
  }
  cat("  \n")
}

BTS_fall_allbox_effic1

BTS_spring_allbox_effic1

Survey weight at age (Atlantis age class)

for(s in names(wtage)){
  cat("  \n##### ",  s,"  \n")
  print(wageplot(wtage[[s]][[1]]))
  cat("  \n")
}

BTS_fall_allbox_effic1

BTS_spring_allbox_effic1

Survey iterpolated weight at age (annual ages)

# annage_wtage read in above 
# annage_wtage <- annage_wtage[[1]] #this still has second list component, diagnostic plot
# annage_wtage <- annage_wtage[[1]] #still in this old version but just removed from function.

for(s in names(annage_wtage)){
  cat("  \n##### ",  s,"  \n")
  print(wageplot(annage_wtage[[s]][[1]][[1]]))
  cat("  \n")
}

BTS_fall_allbox_effic1

BTS_spring_allbox_effic1

Visualize perfectly observed fisheries

We are using the full time series for total catch and weight at age, and a subsample of years 30-35 for compositions. Fishery length and age composition is sampled every output timestep. Fishery effN for length and age sampling was 1000 fish per year (200 per time step), so this is reflected in the composition plots.

Fishery total catch

# compare with true catch, should match
plotC(catchbio_ss, omlist_ss$truecatchbio_ss)

Fishery length composition

lcompsub <- as.data.frame(fish_len_comp_data) %>% filter(time %in% c(150:175)) %>%
  group_by(species) %>%
  group_map(~ plotlen(.x), keep = TRUE)

for(i in 1:length(lcompsub)) {
  print(lcompsub[[i]])
}

Fishery catch at age (Atlantis age class)

acompsub <- as.data.frame(fish_age_comp) %>% filter(time %in% c(150:175)) %>%
  group_by(species) %>%
  left_join(., trueCAAcl) %>%
  #group_map(~ Natageplot(.x), keep = TRUE) # plots only sampled age comp
  group_map(~ Natageplot(.x, effN = 200, truedat = 1), keep = TRUE) # plots merged true age comp


for(i in 1:length(acompsub)) {
  print(acompsub[[i]])
}

Fishery catch at age (annual ages)

acompsub <- as.data.frame(fish_annage_comp) %>% filter(time %in% c(150:175)) %>%
  group_by(species) %>%
  left_join(., trueCAA) %>%
  #group_map(~ Natageplot(.x), keep = TRUE) # plots only sampled age comp
  group_map(~ Natageplot(.x, effN = 200, truedat = 1), keep = TRUE) # plots merged true age comp

for(i in 1:length(acompsub)) {
  print(acompsub[[i]])
}

Fishery weight at age (Atlantis age class)

wageplot(fishwtage)

Fishery interpolated weight at age (annual ages)

#diagnostic plot component still in this old version but just removed from function.
wageplot(fish_annage_wtage[[1]])

Specifying a more realistic survey.

New survey and fishery config files

mssurvey_spring_01.R and mssurvey_fall_01.R configure the fishery independent surveys. We keep spring and fall survey timing, but now we have efficiency defined for each species and a selectivity function for each species. Our surveys exclude offshore areas. These specifications are invented for demonstration purposes, not intended to realistically mimic a particular survey:

# Default survey configuration here has a range of efficiencies and selectivities
# To emulate a range of species in a single multispecies survey
# Also now happens in "spring" and "fall"
# Need to define survey season, area, efficiency, selectivity

# Survey name
survey.name="BTS_spring_nearbox_qmix_selmix"

#Atlantis model timestep corresponding to the true output--now from census_spec.R
timestep <- stepperyr #5

#Which atlantis timestep does the survey run in?--now from census_spec.R
# with 5 output steps per year, 0 is Jan-Feb-midMar, 1 is midMar-Apr-May, 
# 2 is June-July-midAug, 3 is midAug-Sept-Oct, 4 is Nov-Dec (ish)

survey_sample_time <- 1 # spring survey

#The last timestep to sample
total_sample <- noutsteps-1 #495

#Vector of indices of survey times to pull
survey_sample_full <- seq(survey_sample_time,
                          total_sample, by=timestep)

survtime <- survey_sample_full

# survey area
# should return all model areas
# survboxes <- allboxes

# shelf boxes excluding farthest east Barents
survboxes <- c(1:8, 21, 23:31, 33, 39:44, 47:49)

# survey efficiency (q)
# should return a perfectly efficient survey 
# surveffic <- data.frame(species=survspp,
#                     efficiency=rep(1.0,length(survspp)))

# group species to apply different efficiency (all NOBA fish and shark groups)
nontrawl <- c("Sharks_other", "Pelagic_large","Mesop_fish")
pelagics <- c("Pelagic_small","Redfish_other","Mackerel","Haddock",
              "Saithe","Redfish","Blue_whiting","Norwegian_ssh","Capelin")
demersals <- c("Demersals_other","Demersal_large","Flatfish_other","Skates_rays",
               "Green_halibut","North_atl_cod","Polar_cod","Snow_crab")
selflats <- c("Long_rough_dab")

#   define bottom trawl mixed efficiency
ef.nt <- 0.01 # for large pelagics, reef dwellers, others not in trawlable habitat
ef.pl <- 0.2  # for pelagics
ef.dm <- 0.7  # for demersals
ef.fl <- 1.1  # for selected flatfish

# bottom trawl survey efficiency specification by species group
effnontrawl <- data.frame(species=nontrawl, efficiency=rep(ef.nt,length(nontrawl)))
effpelagics <- data.frame(species=pelagics, efficiency=rep(ef.pl,length(pelagics)))
effdemersals <- data.frame(species=demersals, efficiency=rep(ef.dm,length(demersals)))
effselflats <- data.frame(species=selflats, efficiency=rep(ef.fl,length(selflats)))

efficmix <- bind_rows(effnontrawl, effpelagics, effdemersals, effselflats)

surveffic <- efficmix %>%
  filter(species %in% survspp)

# survey selectivity (true age based, flat)
# for annage output uses names(annages) NOT alphabetical survspp
 survselex <- data.frame(species=rep(names(annages), n_annages), #  
                         agecl=unlist(sapply(n_annages,seq)),
                         selex=rep(1.0,sum(n_annages)))

#   mixed selectivity: specify for annual ages 0-10 to 
#   then apply that curve to agecl keeping every nth selectivity index n=NumAgeClassSize
#   ageclsel <- fullsel[seq(NumAgeClassSize, n_annages, length.out=NumCohorts)]
#     flat=1 for large pelagics, reef dwellers, others not in trawlable habitat
#     sigmoid 0 to 1 with 0.5 inflection ~ age 3 for pelagics, reaching 1 at age 5, flat top
#     sigmoid 0 to 1 with 0.5 inflection ~ age 5 for most demersals and flatfish, reaching 1 at age 7, flat top
#     dome shaped 0 to 1 at agecl 6&7 for selected demersals, falling off to 0.7 by agecl 10--didn't do

sigmoid <- function(a,b,x) {
  1 / (1 + exp(-a-b*x))
}

 sp_age <- sp_age %>%
   mutate(n_annages = NumCohorts * NumAgeClassSize) 
# 
# survey selectivity specification for true ages 1-10 by species group--replace for each group in surselex

selnontrawl <- data.frame(species=rep(nontrawl, each=10),
                          agecl=rep(c(1:10),length(nontrawl)),
                          selex=rep(1.0,length(nontrawl)*10))
selpelagics <- data.frame(species=rep(pelagics, each=10),
                          agecl=rep(c(1:10),length(pelagics)),
                          selex=sigmoid(5,1,seq(-10,10,length.out=10)))
seldemersals <- data.frame(species=rep(demersals, each=10),
                           agecl=rep(c(1:10),length(demersals)),
                           selex=sigmoid(1,1,seq(-10,10,length.out=10)))
selselflats <- data.frame(species=rep(selflats, each=10),
                          agecl=rep(c(1:10),length(selflats)),
                          selex=sigmoid(1,1,seq(-10,10,length.out=10)))

selexmix <- bind_rows(selnontrawl, selpelagics, seldemersals, selselflats) 

selexmix <- selexmix %>%
  filter(species %in% survspp) %>%
  rename(selex10 = selex)

survselex <- merge(survselex, selexmix, all = TRUE) %>%
  filter(!is.na(selex)) %>%
  mutate(selex = case_when(!is.na(selex10) ~ selex10,
                           is.na(selex10) ~ selex)) %>%
  select(-selex10)
  
# now apply this for agecl selectivity that matches
# and figure out how to use in wrapper!
#ageclsel <- fullsel[seq(NumAgeClassSize, n_annages, length.out=NumCohorts)]

survselex.agecl <- survselex %>% left_join(sp_age, by=c("species"="Name")) %>%
  group_by(species) %>%
  filter(agecl %in% seq(unique(NumAgeClassSize), 
                        unique(n_annages), 
                        length.out=unique(NumCohorts))) %>%
  mutate(agecl = agecl/unique(NumAgeClassSize)) %>%
  select(species, agecl, selex)

# effective sample size needed for sample_fish
# this is the number of *lengths* per species that are measured on a survey
# effective N of 100000 with no age subsample matches true age comp
surveffN <- data.frame(species=survspp, effN=rep(10000, length(survspp)))

# survey index cv needed for sample_survey_xxx
# cv = 0.1
# surv_cv <- data.frame(species=survspp, cv=rep(0.1,length(survspp)))

# use this constant 0 cv for testing
surv_cv_0 <- data.frame(species=survspp, cv=rep(0.0,length(survspp)))

#   define bottom trawl survey cv by group
cv.nt <- 1.0 # for large pelagics, reef dwellers, others not in trawlable habitat
cv.pl <- 0.4  # for pelagics
cv.dm <- 0.2  # for demersals
cv.fl <- 0.2  # for selected flatfish

# specify cv by species groups
surv_cv_nontrawl <- data.frame(species=nontrawl, cv=rep(cv.nt,length(nontrawl)))
surv_cv_pelagics <- data.frame(species=pelagics, cv=rep(cv.pl,length(pelagics)))
surv_cv_demersals <- data.frame(species=demersals, cv=rep(cv.dm,length(demersals)))
surv_cv_selflats <- data.frame(species=selflats, cv=rep(cv.fl,length(selflats)))

surv_cv_mix <- bind_rows(surv_cv_nontrawl, surv_cv_pelagics, surv_cv_demersals, surv_cv_selflats)

surv_cv <- surv_cv_mix %>%
  filter(species %in% survspp)

# length at age cv for input into calc_age2length function
# function designed to take one cv for all species, need to change to pass it a vector
lenage_cv <- 0.1

# max size bin for length estimation, function defaults to 150 cm if not supplied
maxbin <- 200

# Default survey configuration here has a range of efficiencies and selectivities
# To emulate a range of species in a single multispecies survey
# Also now happens in "spring" and "fall"
# Need to define survey season, area, efficiency, selectivity

# Survey name
survey.name="BTS_fall_nearbox_qmix_selmix"

#Atlantis model timestep corresponding to the true output--now from census_spec.R
timestep <- stepperyr #5

#Which atlantis timestep does the survey run in?--now from census_spec.R
# with 5 output steps per year, 0 is Jan-Feb-midMar, 1 is midMar-Apr-May, 
# 2 is June-July-midAug, 3 is midAug-Sept-Oct, 4 is Nov-Dec (ish)

survey_sample_time <- 3 # fall survey

#The last timestep to sample
total_sample <- noutsteps-1 #495

#Vector of indices of survey times to pull
survey_sample_full <- seq(survey_sample_time,
                          total_sample, by=timestep)

survtime <- survey_sample_full

# survey area
# should return all model areas
#survboxes <- allboxes
# shelf boxes excluding farthest east Barents
survboxes <- c(1:8, 21, 23:31, 33, 39:44, 47:49)

# survey efficiency (q)
# should return a perfectly efficient survey 
# surveffic <- data.frame(species=survspp,
#                     efficiency=rep(1.0,length(survspp)))

# group species to apply different efficiency (all NOBA fish and shark groups)
nontrawl <- c("Sharks_other", "Pelagic_large","Mesop_fish")
pelagics <- c("Pelagic_small","Redfish_other","Mackerel","Haddock",
              "Saithe","Redfish","Blue_whiting","Norwegian_ssh","Capelin")
demersals <- c("Demersals_other","Demersal_large","Flatfish_other","Skates_rays",
               "Green_halibut","North_atl_cod","Polar_cod","Snow_crab")
selflats <- c("Long_rough_dab")

#   define bottom trawl mixed efficiency
ef.nt <- 0.01 # for large pelagics, reef dwellers, others not in trawlable habitat
ef.pl <- 0.1  # for pelagics
ef.dm <- 0.6  # for demersals
ef.fl <- 1.1  # for selected flatfish

# bottom trawl survey efficiency specification by species group
effnontrawl <- data.frame(species=nontrawl, efficiency=rep(ef.nt,length(nontrawl)))
effpelagics <- data.frame(species=pelagics, efficiency=rep(ef.pl,length(pelagics)))
effdemersals <- data.frame(species=demersals, efficiency=rep(ef.dm,length(demersals)))
effselflats <- data.frame(species=selflats, efficiency=rep(ef.fl,length(selflats)))

efficmix <- bind_rows(effnontrawl, effpelagics, effdemersals, effselflats)

surveffic <- efficmix %>%
  filter(species %in% survspp)

# survey selectivity (true age based, flat)
# for annage output uses names(annages) NOT alphabetical survspp
survselex <- data.frame(species=rep(names(annages), n_annages), #  
                        agecl=unlist(sapply(n_annages,seq)),
                        selex=rep(1.0,sum(n_annages)))

#   mixed selectivity: specify for annual ages 0-10 to 
#   then apply that curve to agecl keeping every nth selectivity index n=NumAgeClassSize
#   ageclsel <- fullsel[seq(NumAgeClassSize, n_annages, length.out=NumCohorts)]
#     flat=1 for large pelagics, reef dwellers, others not in trawlable habitat
#     sigmoid 0 to 1 with 0.5 inflection ~ age 3 for pelagics, reaching 1 at age 5, flat top
#     sigmoid 0 to 1 with 0.5 inflection ~ age 5 for most demersals and flatfish, reaching 1 at age 7, flat top
#     dome shaped 0 to 1 at agecl 6&7 for selected demersals, falling off to 0.7 by agecl 10--didn't do

sigmoid <- function(a,b,x) {
  1 / (1 + exp(-a-b*x))
}

sp_age <- sp_age %>%
  mutate(n_annages = NumCohorts * NumAgeClassSize) 
# 
# survey selectivity specification for true ages 1-10 by species group--replace for each group in surselex

selnontrawl <- data.frame(species=rep(nontrawl, each=10),
                          agecl=rep(c(1:10),length(nontrawl)),
                          selex=rep(1.0,length(nontrawl)*10))
selpelagics <- data.frame(species=rep(pelagics, each=10),
                          agecl=rep(c(1:10),length(pelagics)),
                          selex=sigmoid(5,1,seq(-10,10,length.out=10)))
seldemersals <- data.frame(species=rep(demersals, each=10),
                           agecl=rep(c(1:10),length(demersals)),
                           selex=sigmoid(1,1,seq(-10,10,length.out=10)))
selselflats <- data.frame(species=rep(selflats, each=10),
                          agecl=rep(c(1:10),length(selflats)),
                          selex=sigmoid(1,1,seq(-10,10,length.out=10)))

selexmix <- bind_rows(selnontrawl, selpelagics, seldemersals, selselflats) 

selexmix <- selexmix %>%
  filter(species %in% survspp) %>%
  rename(selex10 = selex)

survselex <- merge(survselex, selexmix, all = TRUE) %>%
  filter(!is.na(selex)) %>%
  mutate(selex = case_when(!is.na(selex10) ~ selex10,
                           is.na(selex10) ~ selex)) %>%
  select(-selex10)

# now apply this for agecl selectivity that matches
# and figure out how to use in wrapper!
#ageclsel <- fullsel[seq(NumAgeClassSize, n_annages, length.out=NumCohorts)]

survselex.agecl <- survselex %>% left_join(sp_age, by=c("species"="Name")) %>%
  group_by(species) %>%
  filter(agecl %in% seq(unique(NumAgeClassSize), 
                        unique(n_annages), 
                        length.out=unique(NumCohorts))) %>%
  mutate(agecl = agecl/unique(NumAgeClassSize)) %>%
  select(species, agecl, selex)

# effective sample size needed for sample_fish
# this is the number of *lengths* per species that are measured on a survey
# effective N of 100000 with no age subsample matches true age comp
surveffN <- data.frame(species=survspp, effN=rep(10000, length(survspp)))

# survey index cv needed for sample_survey_xxx
# cv = 0.1
# surv_cv <- data.frame(species=survspp, cv=rep(0.1,length(survspp)))

# use this constant 0 cv for testing
surv_cv_0 <- data.frame(species=survspp, cv=rep(0.0,length(survspp)))

#   define bottom trawl survey cv by group
cv.nt <- 1.0 # for large pelagics, reef dwellers, others not in trawlable habitat
cv.pl <- 0.4  # for pelagics
cv.dm <- 0.2  # for demersals
cv.fl <- 0.2  # for selected flatfish

# specify cv by species groups
surv_cv_nontrawl <- data.frame(species=nontrawl, cv=rep(cv.nt,length(nontrawl)))
surv_cv_pelagics <- data.frame(species=pelagics, cv=rep(cv.pl,length(pelagics)))
surv_cv_demersals <- data.frame(species=demersals, cv=rep(cv.dm,length(demersals)))
surv_cv_selflats <- data.frame(species=selflats, cv=rep(cv.fl,length(selflats)))

surv_cv_mix <- bind_rows(surv_cv_nontrawl, surv_cv_pelagics, surv_cv_demersals, surv_cv_selflats)

surv_cv <- surv_cv_mix %>%
  filter(species %in% survspp)

# length at age cv for input into calc_age2length function
# function designed to take one cv for all species, need to change to pass it a vector
lenage_cv <- 0.1

# max size bin for length estimation, function defaults to 150 cm if not supplied
maxbin <- 200

msfishery_01.R configures the fishery dependent data, which has a higher observation error on the total catch and does not get biological sampling from all areas:

# Default fishery configuration here is a census
fishery.name="nearbox"

# Inherits species from input omlist_ss
fishspp <- omlist_ss$species_ss 

#Number of years of data to pull
nyears <- 50

#Atlantis initialization period in years
burnin <- 30

# fishery output: learned the hard way this can be different from ecosystem outputs
fstepperyr <- if(omlist_ss$runpar$outputstepunit=="days") 365/omlist_ss$runpar$toutfinc


# same time dimensioning parameters as in surveycensus.R
#Vector of indices of catch in numbers to pull (by timestep to sum)
fish_sample_full <- c(0:total_sample)  #total_sample defined in survey configs
fish_burnin <- burnin*fstepperyr+1
fish_nyears <- nyears*fstepperyr
fish_times <- fish_sample_full[fish_burnin:(fish_burnin+fish_nyears-1)]
fish_timesteps <- seq(fish_times[fstepperyr], max(fish_times), by=fstepperyr) #last timestep
#fish_years <- unique(floor(fish_times/fstepperyr)+1) # my original
fish_years <- unique(floor(fish_times/fstepperyr)) #from Christine's new sardine_config.R

fishtime <- fish_times


# fishery sampling area (applies to comps, not total catch)
# should return all model areas, this assumes you see everything that is caught
#fishboxes <- c(0:(omlist_ss$boxpars$nbox - 1))
# same as survey, only getting catch comps from nearby boxes
fishboxes <- c(1:8, 21, 23:31, 33, 39:44, 47:49)

# effective sample size needed for sample_fish
# this effective N is divided by the number of annual timesteps below, so 200 per time
# use as input to the length samples, ages can be a subset
fisheffN <- data.frame(species=survspp, effN=rep(1000, length(survspp)))

#this adjusts for subannual fishery output so original effN is for whole year
fisheffN$effN <- fisheffN$effN/fstepperyr 

# fishery catch cv can be used in sample_fishery_totcatch to
# perfect observation is cv=0
fish_cv <- data.frame(species=survspp, cv=rep(0.1,length(survspp)))

Model boxes show where our survey works now. We are including boxes 1:8, 21, 23:31, 33, 39:44, 47:49 in these surveys, and in fishery comps sampling:

knitr::include_graphics("visualization_Rmd/images/NOBAboxes.png")

The way sampling is set up, effective N for each survey is the number of fish measured for length, age, and average weight at age. We could introduce further realism by using the sample_ages() function to take a subsample for age composition and optionally apply an ageing error matrix. This would require changing the om_comps() wrapper: see atlantisom issue.

New survey and fishery outputs

#NOBAom <- om_init(here("config/NOBA2config.R"))

#NOBAom_ms <- om_species(sppsubset$Name, NOBAom)

NOBAom_ms <- readRDS(file.path(d.name, paste0(scenario.name, "omlist_ss.rds")))

# new wrapper for multiple surveys saves survey output separately
# but returns a list with all surveys, not ideal; make consistent

NOBAom_ms_ind <- om_index(usersurvey = c(here("config/mssurvey_spring_01.R"), 
                                         here("config/mssurvey_fall_01.R")), 
                           userfishery = here("config/msfishery_01.R"),
                           omlist_ss = NOBAom_ms, 
                           n_reps = 1, 
                           save = TRUE)

NOBAom_ms_comp <- om_comps(usersurvey = c(here("config/mssurvey_spring_01.R"), 
                                         here("config/mssurvey_fall_01.R")), 
                           userfishery = here("config/msfishery_01.R"),
                           omlist_ss = NOBAom_ms, 
                           n_reps = 1, 
                           save = TRUE)

NOTE! Need to have selectivities defined both for annual ages and for age classes. Worked in census because flat selectivity was defined for full set of annual ages and it just didn’t use past age 10 for age class output. They were all 1 anyway so it didnt matter. The updated survey config file first specifies selectivity by true age and then backs out an approximate equivalent for age classes. We will use true age classes for the real datasets but wanted to check both outputs here. om_comps() has now changed to take survselex.agecl as input to the age class comps and survselex as input to the full comps.

Both om_index() and om_comps() wrappers now name fishery output based on a fishery.name in the fishery config file.